Data Preprocessing

Load and annotate data

# Load csvs
datExpr = read.csv('./../Data/Gandal/RNAseq_ASD_datExpr.csv', row.names=1)
datMeta = read.csv('./../Data/Gandal/RNAseq_ASD_datMeta.csv')

# Group brain regions by lobes
datMeta$Brain_Region = as.factor(datMeta$Region)
datMeta$Brain_lobe = 'Occipital'
datMeta$Brain_lobe[datMeta$Brain_Region %in% c('BA4_6', 'BA9', 'BA24', 'BA44_45')] = 'Frontal'
datMeta$Brain_lobe[datMeta$Brain_Region %in% c('BA3_1_2_5', 'BA7')] = 'Parietal'
datMeta$Brain_lobe[datMeta$Brain_Region %in% c('BA38', 'BA39_40', 'BA20_37', 'BA41_42_22')] = 'Temporal'
datMeta$Brain_lobe=factor(datMeta$Brain_lobe, levels=c('Frontal', 'Temporal', 'Parietal', 'Occipital'))

# Remove '/' from Batch variable: (It is recommended (but not required) to use only letters, numbers, 
# and delimiters '_' or '.', in levels of factors as these are safe characters for column names in R
datMeta$Batch = gsub('/', '.', datMeta$RNAExtractionBatch) %>% as.factor

# Transform Diagnosis into a factor variable
datMeta$Diagnosis_ = factor(datMeta$Diagnosis_, levels=c('CTL','ASD'))

Check sample composition

Data description taken from the dataset’s synapse entry: RNAseq data was generated from 88 postmortem cortex brain samples from subjects with ASD (53 samples from 24 subjects) and non-psychiatric controls (35 samples from 17 subjects), across four cortical regions encompassing all major cortical lobes – frontal, temporal, parietal, and occipital. Brain samples were obtained from the Harvard Brain Bank as part of the Autism Tissue Project (ATP).

print(paste0('Dataset includes ', nrow(datExpr), ' probes from ', ncol(datExpr), ' samples belonging to ', length(unique(datMeta$Subject_ID)), ' different subjects.'))

## [1] "Dataset includes 63682 probes from 88 samples belonging to 41 different subjects."

Diagnosis distribution: There are more ASD samples than controls

table(datMeta$Diagnosis_)

## 
## CTL ASD 
##  35  53

Brain region distribution: All regions seem to be balanced

table(datMeta$Brain_lobe)

## 
##   Frontal  Temporal  Parietal Occipital 
##        22        20        23        23

Diagnosis and brain region seem to be balanced except for the frontal lobe, where there are more control samples than ASD ones.

table(datMeta$Diagnosis_, datMeta$Brain_lobe)

##      
##       Frontal Temporal Parietal Occipital
##   CTL      13        6        8         8
##   ASD       9       14       15        15

Sex distribution: There are many more Male samples than Female ones

table(datMeta$Sex)

## 
##  F  M 
## 15 73

Age distribution: Subjects between 5 and 60 years old with a mean close to 30

summary(datMeta$Age)

##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##    2.00   17.00   28.00   29.74   41.75   60.00

Annotate probes with BioMart information

Processing group

All the information we have is the Brain Bank, and although all the samples were obtained from the Autism Tissue Project, we don’t have any more specific information about who preprocessed each sample

table(datMeta$Brain_Bank)

## 
## ATP 
##  80

Date of processing

There are two different dates when the data was procesed

table(datMeta$RNAExtractionBatch)

## 
## 10/10/2014  6/20/2014 
##         53         27

Luckily, there doesn’t seem to be a correlation between the batch surrogate and the objective variable, so the batch effect will not get confused with the Diagnosis effect

table(datMeta$RNAExtractionBatch, datMeta$Diagnosis_)

##             
##              CTL ASD
##   10/10/2014  24  29
##   6/20/2014   11  16

*All the samples from each subject were processed on the same day (makes sense, otherwise they wound need to freeze the samples)

Samples don’t seem to cluster together that strongly for each batch, although there does seem to be some kind of relation

h_clusts = datExpr %>% t %>% dist %>% hclust %>% as.dendrogram

create_viridis_dict = function(){
  min_age = datMeta$Age %>% min
  max_age = datMeta$Age %>% max
  viridis_age_cols = viridis(max_age - min_age + 1)
  names(viridis_age_cols) = seq(min_age, max_age)
  
  return(viridis_age_cols)
}
viridis_age_cols = create_viridis_dict()

dend_meta = datMeta[match(substring(labels(h_clusts),2), datMeta$Dissected_Sample_ID),] %>% 
            mutate('Batch' = ifelse(RNAExtractionBatch=='10/10/2014', '#F8766D', '#00BFC4'),
                   'Diagnosis' = ifelse(Diagnosis_=='CTL','#008080','#86b300'), # Blue control, Green ASD
                   'Sex' = ifelse(Sex=='F','#ff6666','#008ae6'),                # Pink Female, Blue Male
                   'Region' = case_when(Brain_lobe=='Frontal'~'#F8766D',        # ggplot defaults for 4 colours
                                        Brain_lobe=='Temporal'~'#7CAE00',
                                        Brain_lobe=='Parietal'~'#00BFC4',
                                        Brain_lobe=='Occipital'~'#C77CFF'),
                   'Age' = viridis_age_cols[as.character(Age)]) %>%             # Purple: young, Yellow: old
            dplyr::select(Age, Region, Sex, Diagnosis, Batch)
h_clusts %>% set('labels', rep('', nrow(datMeta))) %>% set('branches_k_color', k=9) %>% plot
colored_bars(colors=dend_meta)

There seems to be a different behaviour by batch mainly in the first and third principal components

pca = datExpr %>% t %>% prcomp
summary(pca)$importance[,1:3]

##                                 PC1          PC2         PC3
## Standard deviation     273215.31037 171051.12817 98195.28432
## Proportion of Variance      0.57415      0.22504     0.07416
## Cumulative Proportion       0.57415      0.79920     0.87336

plot_data = data.frame('ID'=colnames(datExpr), 'PC1' = pca$x[,1], 'PC2' = pca$x[,2], 'PC3' = pca$x[,3]) %>% 
            mutate('ID'=substring(ID,2)) %>% left_join(datMeta, by=c('ID'='Dissected_Sample_ID')) %>% 
            mutate('Batch'=RNAExtractionBatch) %>% dplyr::select('PC1','PC2','PC3','Batch')

plot_data %>% ggpairs(progress=FALSE, aes(colour=Batch, fill=Batch, alpha=0.3)) + theme_minimal()

## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

rm(pca, plot_data)

Comparing the mean expression of each sample by batch we can see the batch effect is reflected here as well

plot_data_b1 = data.frame('Mean'=colMeans(datExpr[,datMeta$RNAExtractionBatch=='10/10/2014']), 'Batch'='10/10/2014')
plot_data_b2 = data.frame('Mean'=colMeans(datExpr[,datMeta$RNAExtractionBatch=='6/20/2014']), 'Batch'='6/20/2014')

plot_data = rbind(plot_data_b1, plot_data_b2)
mu = plot_data %>% group_by(Batch) %>% dplyr::summarise(BatchMean=mean(Mean))

ggplotly(plot_data %>% ggplot(aes(x=Mean, color=Batch, fill=Batch)) + geom_density(alpha=0.3) + 
         geom_vline(data=mu, aes(xintercept=BatchMean, color=Batch), linetype='dashed') +
         ggtitle('Mean expression by sample grouped by Batch') + scale_x_log10() + theme_minimal())

rm(plot_data_b1, plot_data_b2, plot_data)

In conclusion: Date of extraction works as a surrogate for batch effect

Save filtered and annotated dataset

save(datExpr, datMeta, datProbes, file='./../Data/Gandal/filtered_raw_data.RData')

And repeating Normalisation

counts = datExpr %>% as.matrix
rowRanges = GRanges(datProbes$chromosome_name,
                  IRanges(datProbes$start_position, width=datProbes$length),
                  strand=datProbes$strand,
                  feature_id=datProbes$ensembl_gene_id)
se = SummarizedExperiment(assays=SimpleList(counts=counts), rowRanges=rowRanges, colData=datMeta)
dds = DESeqDataSet(se, design = ~ Batch + Diagnosis_) # Including batch effect

# Perform DEA
dds = DESeq(dds)

## estimating size factors

## estimating dispersions

## gene-wise dispersion estimates

## mean-dispersion relationship

## final dispersion estimates

## fitting model and testing

## -- replacing outliers and refitting for 145 genes
## -- DESeq argument 'minReplicatesForReplace' = 7 
## -- original counts are preserved in counts(dds)

## estimating dispersions

## fitting model and testing

DE_info = results(dds, lfcThreshold=log2(1.2), altHypothesis='greaterAbs')

# Perform vst
vsd = vst(dds)

datExpr_vst = assay(vsd)
datMeta_vst = colData(vsd)
datProbes_vst = rowRanges(vsd)

rm(counts, rowRanges, se, vsd)

This plot remains stable

meanSdPlot(datExpr_vst, plot=FALSE)$gg + theme_minimal()

This one looks better now. The valley found in the original density plot of the mean expression of the probes seems to still be present here beecause there seem to be two clouds main of points

plot_data = data.frame('ID'=rownames(datExpr_vst), 'Mean'=rowMeans(datExpr_vst), 'SD'=apply(datExpr_vst,1,sd))

plot_data %>% ggplot(aes(Mean, SD)) + geom_point(color='#0099cc', alpha=0.05) + 
              scale_x_log10() + scale_y_log10() + theme_minimal()

rm(plot_data)

datExpr = datExpr_vst
datMeta = datMeta_vst %>% data.frame
datProbes = datProbes_vst

print(paste0('After filtering, the dataset consists of ', nrow(datExpr), ' probes and ', ncol(datExpr), ' samples'))

## [1] "After filtering, the dataset consists of 30112 probes and 80 samples"

rm(datExpr_vst, datMeta_vst, datProbes_vst)

Performing Batch Correction using ComBat

datExpr = datExpr %>% as.matrix %>% ComBat(batch=datMeta$Batch)

## Found2batches

## Adjusting for0covariate(s) or covariate level(s)

## Standardizing Data across genes

## Fitting L/S model and finding priors

## Finding parametric adjustments

## Adjusting the Data

Plots changed although I’m not sure if it’s for better or for worse?

pca = datExpr %>% t %>% prcomp
summary(pca)$importance[,1:3]

##                             PC1      PC2      PC3
## Standard deviation     35.78914 29.62096 21.34238
## Proportion of Variance  0.19990  0.13694  0.07109
## Cumulative Proportion   0.19990  0.33684  0.40793

plot_data = data.frame('ID'=colnames(datExpr), 'PC1' = pca$x[,1], 'PC2' = pca$x[,2], 'PC3' = pca$x[,3]) %>% 
            mutate('ID'=substring(ID,2)) %>% left_join(datMeta, by=c('ID'='Dissected_Sample_ID')) %>% 
            mutate('Batch'=RNAExtractionBatch) %>% dplyr::select('PC1','PC2','PC3','Batch')

plot_data %>% ggpairs(progress=FALSE, aes(colour=Batch, fill=Batch, alpha=0.3)) + theme_minimal()

## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

rm(pca, plot_data)

Both batches now have almost the same mean expression

plot_data_b1 = data.frame('Mean'=colMeans(datExpr[,datMeta$RNAExtractionBatch=='10/10/2014']), 'Batch'='10/10/2014')
plot_data_b2 = data.frame('Mean'=colMeans(datExpr[,datMeta$RNAExtractionBatch=='6/20/2014']), 'Batch'='6/20/2014')

plot_data = rbind(plot_data_b1, plot_data_b2)
mu = plot_data %>% group_by(Batch) %>% dplyr::summarise(BatchMean=mean(Mean))

ggplotly(plot_data %>% ggplot(aes(x=Mean, color=Batch, fill=Batch)) + geom_density(alpha=0.3) + 
         geom_vline(data=mu, aes(xintercept=BatchMean, color=Batch), linetype='dashed') +
         ggtitle('Mean expression by sample grouped by Batch') + scale_x_log10() + theme_minimal())

rm(plot_data_b1, plot_data_b2, plot_data)

Save preprocessed dataset

save(datExpr, datMeta, datProbes, DE_info, dds, file='./../Data/Gandal/preprocessed_data.RData')

Load SFARI Genes and annotate them with ensembl IDs

# Load SFARI information
SFARI_genes = read_csv('./../Data/SFARI/SFARI_genes_01-15-2019.csv')

# Get ensemble IDS for SFARI genes
mart = useMart(biomart='ENSEMBL_MART_ENSEMBL',
               dataset='hsapiens_gene_ensembl',
               host='feb2014.archive.ensembl.org') ## Gencode v19

gene_names = getBM(attributes=c('ensembl_gene_id', 'hgnc_symbol'), filters=c('hgnc_symbol'), 
                   values=SFARI_genes$`gene-symbol`, mart=mart) %>% 
                   mutate('gene-symbol'=hgnc_symbol, 'ID'=as.character(ensembl_gene_id)) %>% 
                   dplyr::select('ID', 'gene-symbol')

SFARI_genes = left_join(SFARI_genes, gene_names, by='gene-symbol') %>% 
              distinct(ID, .keep_all = TRUE)

write.csv(SFARI_genes, './../Data/SFARI/SFARI_genes_with_ensembl_IDs.csv', row.names=F)